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We describe in detail a recently proposed lattice-Boltzmann model for simulating 
flows with multiple phases and components. In particular, the focus is on the 



o 

^ modeling of one-component fluid systems which obey non-ideal gas equations of 

^ state and can undergo a liquid-gas type phase transition. The model is shown to be 

I momentum-conserving. From the microscopic mechanical stability condition, the 

C densities in bulk liquid and gas phases are obtained as functions of a temperature-like 

o 

^ parameter. Comparisons with the thermodynamic theory of phase transition show 

that the LBE model can be made to correspond exactly to an isothermal process. 
The density profile in the liquid-gas interface is also obtained as function of the 
temperature-like parameter and is shown to be isotropic. The surface tension, which 
can be changed independently, is calculated. The analytical conclusions are verified 
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by numerical simulations. 
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I. INTRODUCTION 

Since the Lattice Gas Automaton (LGA) was introduced [1-3] as an alternative method for 
simulating fluid flows, a great amount of effort has been devoted to developing models for various 
physical systems [4]. In a LGA model, both space and time are discrete and particles move along 
links on a regular lattice and collide with each other on each lattice site according to some properly 
designed rules and conservation laws. A set of fluid equations can be derived describing the 
macro-dynamics of the system which, under proper limitations, lead to incompressible Navier- 
Stokes equations. 

Although the LGA offers a simple and efficient tool for simulating some fluid flows, it is difficult 
to achieve high accuracy with the LGA because of its statistical noise. Moreover, it usually suffers 
from other problems including a lack of Galilean invariance and velocity-dependent pressure, 
unless some special procedures are taken to eliminate them [5]. To overcome these problems, an 
alternative approach, known as the lattice Boltzmann Equation (LBE) method [6-8], was derived 
from the LGA. In this method, rather than following dynamics of the particles, a Boltzmann 
equation governing the evolution of the particle distribution function is directly solved. The 
solution, after some simple mapping, can be shown to obey the desired fluid equations under 
certain conditions [9,10]. While having inherited the major advantages of the LGA in simulating 
fluid flows, such as the computational simplicity and the ease in handling complex boundary 
conditions, the LBE method is free of the problems mentioned above. 

There have been LGA models proposed for simulating flows with multiple thermodynamic 
phases. One of the first models was constructed by Chen et al. [11] in which the particles interact 
according to a nearest-neighbor potential. The thermodynamic properties are equivalent to that of 
the Ising model so that phase transition and the critical point is exactly known [12] . However, the 
momentum exchange in this model does not correspond to that in a real fluid. It can not be used 
for simulating realistic fluid flows. 

Alternatively, Appert et al. [13-15] proposed a momentum conserving LGA model for simu- 
lating the liquid-gas phase transition. In their model, particles over several lattice sites can interact 



with each other. As the range of the interaction exceeds a certain critical value, a liquid-gas type of 
phase transition is observed. Nevertheless, this model suffers from the shortcomings of the usual 
LGA models, especially the lack of Galilean invariance which becomes severe when there is a 
large density variation. Moreover, the surface tension in this model depends on the orientation of 
the interface with respect to the lattice structure [15]. 

In a previous paper [16], we presented the basic idea of a new LBE model for modeling flows 
with multiple thermodynamic phases based on nearest- neighbor interactions. In addition to the 
advantages of the LBE models, it offers some features highly desirable in simulating multi-phase 
flows. First of all, it is computationally efficient since for most purposes, the non-locality can be 
restricted to nearest- neighbor interactions. The additional computation involved in such a non-local 
interaction is equivalent to that of a streaming step, which can be efficiently implemented on parallel 
computers. Secondly, it can handle fluid systems composed of arbitrary number of components, 
each of which can have different mass and transport coefficients. Third, the equation of state is 
exactly expressed in terms of the inter-particle potential, and can be tuned to precisely match any 
given functional form. When the equation of state is properly chosen, a liquid-gas phase transition 
can occur. Critical phenomena can be studied since the critical point is analytically calculable. 

In this paper, we focus our attention to the application of this model to one-component systems 
and discuss in detail the equilibrium properties when a liquid-gas phase transition is allowed. 
In Section II, we present the details of the inter-particle interactions and show the momentum 
conservation in a global sense. The liquid-gas phase transition is shown to occur when the inter- 
particle potential is properly chosen and a temperature-like parameter is below a critical value. 
In Section III, the two-phase coexistence curve is calculated analytically from the microscopic 
mechanical balance condition. The pressure tensor is given in analytical form, from which the 
density profile across the liquid-gas interface is obtained. The surface tension is also calculated as 
a function of the temperature-like parameter. Comparison with thermodynamic predictions reveals 
that this LBE model corresponds to an isothermal PVT system when the inter-particle potential is 
chosen to be in a particular form. In Section IV, the analytical results are verified by numerical 
computation. The agreement between the theory and the simulation is found to be excellent. The 



density profile is also shown to be isotropic with respect to the underlying lattice structure. Finally 
in Section V, we offer further discussion and suggest some future research directions. 

II. LATTICE-BOLTZMANN MODEL WITH PAIRWISE INTER-PARTICLE POTENTIALS 

We recall the one-component lattice Boltzmann equation with the BGK [9,17] collision term: 

na(x + e^, t + 1) - na(x, t) = [^^(x, t) - <^(x, t)] , (1) 

T 

where the distribution function ^^(x, t) at lattice site x denotes the particle population moving in 
the direction e^; n^{x, t) is a prescribed equilibrium distribution function; and r is the collision 
time [9]. It has been shown [9,10] that if the equilibrium distribution function is chosen to be 



<^(x) = n(x) [i^ + ^e„ • u + ^^e^e, : uu - ^u^] , a=l,---,b; 



nn^fxl = n(x.) \do - 4u^ 



(2) 



at long wave length limit, Navier-Stokes equations with an ideal gas equation of state 

P = ^ n (3) 

will be obtained from the kinetic equation (1). The corresponding kinematic shear viscosity, 
u, is c^{t — \)/{D + 2) and the bulk viscosity, rj, is c?oc^(t — \)/D. In the equations above, 
'^(x) = J2a'>T'a{^) and u = J2a^ana{^)/n{x.) are the particle number and the fluid velocity at 
lattice site x. b is the number of total links connecting a lattice site to its nearest neighbors. 
The velocity vectors are represented by e^ (a = 1, • • • , 6) with magnitudes equal to c, the lattice 
constant divided by the time step. D is the dimension of the chosen lattice [3]. The constant, 
do, is the equilibrium fraction of particles with zero speed. It can be easily verified that with the 
equilibrium distribution function, nl^, so chosen, the BGK collision operator on the right hand side 
of Eq. (1) conserves the total particle number, n(x), and the total momentum, X^„ eana(x), at each 
site. 

To simulate non-ideal gases and their mixtures, long-range interactions between particles 
must be included. Since particles then exchange momentum through the long-range attractive or 



repulsive forces in addition to short range collisions, site-wise momentum conservation must be 
abandoned. In a previous paper [16], we suggested the inclusion of effects of the inter-particle 

forces in a system composed of S components by defining a potential of the following form, 

s s 
F(x,x') = E E G'.^(x,x')V''^K(x)]V^K(x')], (4) 

o-=1ct=1 

where, n'^ and n'^ are the number density of components a and a. Since in a lattice model, particles 
reside on lattice sites with fixed distance between each other, the number density on each site also 
determines the average distance between particles. The potential is then made to be proportional 
to the products of local "effective masses" ip'^{n'^), which are functions of local densities and are 
defined individually for different components. Although the exact mapping is currently unknown, 
the forms of the ip's control the detailed nature of the interaction potential and determine the 
equation of state of the system. The forces between different components can be either attractive 
or repulsive and can have different strengths, as defined by the Green's function, G^a (x, x'), which 
is reduced to (jo-ct(x — x') in homogeneous systems. 

In the case of a one-component system with nearest-neighbor interaction, the Green's function 
is further reduced to a single number Q, namely 



^(x - x') = < 



Ix - x'l > c 

(5) 

/^ I /I 

y X — x' = c , 



which measures the strength of the interactions among the particles on the nearest neighboring 
sites. The long-range force due to the potential Eq. (4) acting on the particles at site x is then 

b 
F(x,t) = -^V(x,t)EV'(x + e„,t)e„, (6) 

a=l 

where ■0(x, t) = ■0['^(x, t)]. To reflect the momentum change caused by this force at each time 
step, we let u in the equilibrium distribution function given by Eq. (2) be 

1 



u 



n X 



b b 

J2 ea?T-a(x) - Tgip{-K, t) ^ ■0(x + 6^, t)ea 
.a=l a=l 



(7) 



The second term is the momentum change due to the long-range force between particles at site x 
and its nearest neighbors. 



With this choice of interaction, although the momentum is not conserved by the BGK collision 

operator at each site, it can be easily shown that the total momentum is indeed conserved. At each 

step, summing over all the lattice sites, we have the following net momentum change: 

b 
AP = -gJ2J2 V'(x)V'(x + e„)e„. (8) 

X a=l 

By replacing the dummy variable e^ with — e^, the above is also equal to, 

AP = ^EEV'(x)V'(x-e„)e„. (9) 

X a=l 

If the boundaries do not introduce any momentum flux into the system, then x can be viewed as a 

dummy variable. Let x' = x — e^ and then drop the primes, the above equation becomes 

b 
AP = ^ E E V'(x + e„)V'(x)e„ = - AP, (10) 

X a=l 

and hence, AP = 0. In another words, no net momentum change is incurred by the interaction 
between particles introduced above. Same conclusion can be proved for systems with more than 
one component. 

With the new equilibrium distribution function given by Eqs. (2) and (7), following the same 
Chapman-Enskog procedure [3], we have previously shown [16] that the Navier-Stokes equation 
with the equation of state 

b 



c2 



■l-do)n + -g',p\n] 



(11) 



can be obtained. Fluids with a quite general class of equations of state can then be modeled by 
selecting the function '4){n). Once the function, ■0, is properly chosen, this equation of state will 
exhibit many essential features of a liquid-gas phase transition. In the limit of weak non-local 
force, namely Q ^ 0, Eq. (11) approaches the equation of state of an ideal-gas given by Eq. (3) 
with a constant temperature proportional to ^(1 — c?o). As the parameter —(1 — c?o)/^ decreases 
below a critical value, (either by increasing the inter-particle force or increasing the rest particles 
at each site), pressure p becomes a non-monotonic function of n, qualitatively similar to the van 
derWaals equation of state for many choices of ■!/;, e.g., ■0(n) = 1 — exp(— n) [16]. The unphysical 



negative compressibility dp/dn corresponds to a thermodynamic instability and the system will 

segregate into a denser (liquid) phase and a lighter (gas) phase. The critical value of (1 — do)/Q 

will be given by the following two equations: 

dp c^g (I- do A d^p hc^g ( „ „x 

^ = -^ — r^ + HV =0, and -4 = — ^ U^" + V' ) = 0. (12) 



dn D \ g ' ^^ J ' 5^2 J) 

The equation of state (11) gives constant relations between the pressure and the density, similar 
to the isotherm in a PVT system. The parameter —(1 — c?o)/^ formally plays the same role in 
this LBE model as the temperature in the van der Waals theory of phase transition. However, it 
should be noted that the original kinetic equations (1), (2) and (7) do not imply explicit energy 
conservation. The existence of an energy-like quantity is unknown at this point. As a consequence, 
this LBE model does not have a well defined statistical mechanical temperature, and in general, 
the correspondence between the LBE model and an isothermal process cannot be established. A 
more detailed discussion will be given later in Section III. 

III. EQUILIBRIUM PROPERTIES 

It is of great interest to calculate the coexistence curve for a two-phase system. For the PVT 
system, the coexistence curve is given by a procedure known as the Maxwell construction, which 
requires that the Gibbs potential is minimized around the phase transition where it is triple-valued. 
But for the current lattice Boltzmann model, there is not a priori an explicit energy conservation 
relation. Whether one can find an effective free energy for this system is still unknown. However, 
the simplicity of the form of the inter-particle force in our model enables us to obtain, from the 
microscopic mechanical balance conditions, the coexistence curve as well as the density profile 
across the liquid-gas interface. The surface tension can also be calculated once the pressure tensor 
and the density profile are known. 

We consider the one-component lattice Boltzmann model with nearest neighbor interaction as 
described in last section. The momentum flux tensor comprises a kinetic term due to the free 
streaming of the particles and an additional potential term from the long-range inter-particle force. 
For the inter-particle force given by Eq. (6), the momentum flux tensor is 

8 



n = Xl^-aeaBa + -■0(x) ^VCX + ea)eaea, (13) 

a a 

where the first term is the conventional kinetic term and the second is the potential term. We 
consider the system in equilibrium, namely ^^(x, t) = ^^(x). We further assume there is no net 
mass transfer along any of the links connecting two lattice sites. This can be expressed as 

na{x.) = na{^ + ea), for e^ = -e^. (14) 

For a homogeneous density distribution, the fluid velocities before and after the collision step equal 
to each other and vanish when the system is in equilibrium. But they are not equal in the interface 
region due to the momentum change included in the collision step. Neither velocity is zero in this 
interface zone. However, the magnitude of the fluid velocity of the equilibrium distribution given 
by Eqs. (2) and (7) can be calculated. Multiply the lattice Boltzmann equation (1) by e^ and sum 
over a. Taking into consider Eq. (14), we have, for steady state, the following equation at each site 

—nu — nu = (nu — nu^^) , (15) 

T 

where u^^ = J2a ^an^ is the fluid velocity of the equilibrium distribution and 

nu^^ = nu — C/T'0 ^ ea'0(x + e^) ~ nu — Qt——iJj\/iJj. (16) 

We can then solve for u^^ from the equations above and obtain the following equation, 

/ 1\ c^hg 

Substituting the equilibrium distribution given by Eqs. (2) and (17) into Eq. (13), and using the 
following Taylor expansion of ■0(x + e^), 

V'(x + e„) = V'(x) + VV' • e„ + - VVV' : e„e„ + • • • , (18) 



we obtain the pressure tensor 

p = -c^n + ^V + ^ — rV'VV I + 

^ [ D 2D ^ AD{D + 2) ^ ^J 

55P^^ VV^ + (r - 5) ^^^V^V^. (19) 



where I is the unit tensor. In a homogeneous fluid where all the derivatives vanish, the pressure 
tensor is isotropic, and once again we obtain the equation of state Eq. (11). 

In order to calculate the density profile across a liquid-gas interface, we consider a flat interface 
coincident with the x-y plane separating the two phases. The densities in the bulk liquid and gas 
phases are ni and Ug respectively. We choose the point where n = (n/ + ng)/2 as the origin of the 
z-axis without lose of generality. All the variables vary in the z direction only, so that n = n{z) 
and d/dx = d/dy = 0. In this situation, the pressure tensor p becomes diagonal: 

p(x) = Pxx('Z)exex + Pyy{z)eyey + Pzz{z)ezez, (20) 

where its transverse components pxx{z) = Pyy{z) (= Pt{z)) because of the symmetry. The normal 
component pn{z) = Pzz{z) = Po must be a constant following the mechanical balance condition 
V • p = 0. The constant, pQ, is the hydrodynamic pressure in either bulk phase. For a low-viscosity 
fluid, we can ignore the last term in Eq. (19), since r — ^ ~ i/ <C 1. The normal component of the 
pressure tensor, p^z , can then be further reduced to 



I — do 2 c^bg 2 ^hc'Q 



2 



r{T\ +V.'' 



n 



(21) 



dz J dz^ 

where'0' = dtp /dn and tp" = d^tp/dn^. In either of the liquid and gas phases far from the interface, 
the pressure po satisfies the following relation 

Po = -^—c ni + ^^V (^0 = -^-c Ug + ^^V {ng). (22) 

This is simply the equation of state written for the two bulk phases. Eq. (21), coupled with Eq. (22) 
through the unknown constant po, is the equation from which the density profile n[z)h to be solved 
with the boundary condition dnidz = at 2; = ±00. A simple change of variable can reduce it 
further to obtain a formal solution. Let [dnj dzY = y, and notice ^ = ^^, then Eq. (21) will be 
transformed to the following first-order differential equation of y: 

l-do 2 c^bg 2 3c^6g ip d 2^ /o'5^ 

where '4){n) is a known function of n and all the derivatives are with respect to n. By direct 
integration, we obtain the following solution: 
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SDjD + 2) ff I -do 2 cHg ,\ # 



4(i) + 2) 



,2,/,/2 



3c^ip 



1 (^" " 


1 - C?o 2 

- ^ cn 


2i}po 


2(1 -do) 



\ V^ 2(1 -do) f 
n\ Imp - -^ H 7^^ / Inipdn 



c^bg hg j "^ 1 bg 

where the pressure po and the integration constant will be determined by Eq. (22) and the boundary 
conditions are y{ni) = y{ng) = 0. To be compatible with these boundary conditions, we must 
have 



Jng \ 



PO - -^-c n - ^^V 1 :^dn = 0. (25) 



From Eqs. (25) and (22), we are able to determine the pressure po and the densities ni and rig in 
the bulk liquid and gas phases for any given parameter (1 — do)/g. This defines the coexistence 
curve of this LBE model. 

Once the pressure, po, and the densities, rig and ni, are determined, y{n) = {dn/dzY is 
completely known. The function y{n), and the densities rig and ni, depend only upon a single 
parameter (1 — do)/g. The density profile n{z) across the liquid-gas interface is then given by the 
following equation: 

r y{n)-'^^dn = z, (26) 

J{Tlg+Tll)/2 

which in general can not be expressed in terms of elementary functions except for some very 
special choices of ■0(n). 

The calculation of the surface tension is straight forward. By definition [18], in the case of a 
flat interface, the surface tension can be calculated from the components of the pressure tensor as 
the following integral across the interface 

/oo 
(po -pT)dz. (27) 

-oo 

Using the expressions for the pressure tensor given by Eq. (13), we have 

cHg f°° ,d^ip 



a 



r°° a w 

/ '^^dz. (28) 

J-oo dz^ 



2D{D + 2) 

Integrating by parts and considering the boundary condition that dtp/dz = Oatz = ±00, we obtain 
the following equation 

11 



a c% 



■ / i;%{n)Y^^dn, (29) 



g 1D{D + 2) 

where the right hand side is a function of (1 — d())IQ only. Since n{z) also depends only on 
(1 — do)/G, we are able to change the surface tension independently of the density profile. 

It is interesting to compare the above results with the classic thermodynamic theory of phase 
transition in PVT systems [19,20]. If we take Eq. (11) as an isotherm of a PVT system, the 
Maxwell equal-area construction (which requires /J|^(po — p)dv = 0, where f ~ 1/n is the molar 
volume and vi and Vg are its values in the liquid and gas phases) yields the following equation for 
the coexistence curve. 



f'i 



«-^^^"-w^'-)>=«- <^°> 



Compared with Eq. (25), we find that the coexistence curve of the LBE model will not agree with 
the thermodynamic theory unless we chose the function ■0 to have the form 

Tp{n) = ■!/'oexp(-no/n), (31) 

where ■0o and no are arbitrary constants. This is not surprising since there is not an energy- 
conservation relation guaranteed by the kinetic equation and a well-defined temperature in this 
LBE model. Therefore the LBE model does not necessarily correspond to an exact isothermal 
process. Nevertheless, if we indeed choose ■0 to have the form of Eq. (31), the behavior of the 
LBE model will be consistent with that of an isothermal process, and the parameter ( 1 — do)/Q can 
be used as a temperature scale. We can also derive from the "isotherm" (11) an effective Gibbs 
potential which is always minimum near phase transition, even though we do not have energy 
conservation in the model. 

IV. SIMULATION VERIFICATION 

In this section, we verify the results of Sec. Ill by solving Eqs. (22), (24), (25), (29) numerically 
and comparing the solutions with simulation results. We chose the function ■0(n) to be 1 — exp(— n), 
the same as in the previous paper [16]. On a two-dimensional hexagonal lattice, the critical point 
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of the system can then be calculated from Eq. (12) as —(1 — c?o)/^ = 1.5, at the critical density 
n = ln2~ 0.6931. 

Numerically solving the density profiles from the equations derived in last section is a com- 
plicated procedure. Firstly for a given value of —(1 — do)/Q, the integral on the left hand side of 
Eq. (25) is a function oipo, where the integration limits Ug and ni have to be solved from Eq. (22) 
for given po- We can then solve po together with Ug and ni by finding the root of this function. 
The NAG routine, C05AGE, was utilized to find the root of an arbitrarily given function. Once rig, 
ni and po are obtained, y{n) was found by integrating Eq. (24) using the IMSL routine, QDAGS. 
The density profile, n{z), is then easily found using {dn/dzY = y and employing IMSL routine, 
IVPBS. Integrating Eq. (29), we obtain the surface tension. The whole calculation takes a few 
seconds on a Cray-YMP C90. 

The simulations were carried out on a two-dimensional hexagonal lattice. In most of the 
calculations, the collision time, r, is chosen to be 0.6 unless otherwise specified. We first verify 
the numerically calculated coexistence curve and the density profiles by simulations on a 64 x 256 
lattice. Periodic boundary conditions are used in both directions. The initial density distribution 
was so set up that the density in half of the domain is higher than that in the other half. When the 
system reaches equilibrium after about 10^ iterations, a flat interface parallel to the shorter edge of 
the rectangular domain forms in the middle. The densities in the bulk phases are then measured 
for different values of —(1 — do)/Q and plotted in Fig. 1. The numerical solution of Eqs. (25), 
shown as the solid line, is found to be in very good agreement with simulation results. In Fig. 1, we 
also demonstrate the deviation from the coexistence curve predicted by the Maxwell construction, 
shown as the dashed line. This curve is obtained by solving Eq. (30) instead of Eq. (25) with the 
same procedure. The difference is small but noticeable, especially when —(1 — do)/Q goes far 
below the critical value. Of cause, this difference will not exist if the effective mass '4){n) is chosen 
to have the form given in Eq. (31). 

The density profiles for several different values of the temperature-like parameter below its 
critical value are displayed in Fig. 2. To examine the isotropy of the surface tension, for each value 
of — (1 — c?o)/^, we performed two simulations in which the angles between the normal vector of 

13 



the interface and the lattice links are 30° and 0° respectively [15]. The measured density profiles 
in both cases are displayed. No significant difference between these two cases is found. The 
isotropy of surface tension in this model is also evident as the previously displayed bubbles [16] 
are circular. Here, all the measurements were made at one cross section and at one time without 
any kind of average. The lines shown in Fig. 2 are the density profiles solved from Eqs. (21) and 
(22). The analytical predictions are again in excellent agreement with the simulation results. It 
is also seen that when the temperature-like parameter —(1 — do)IQ approaches the critical point, 
the distinctions between the two phases become smaller and will eventually vanish at the critical 
temperature. 

Plotted in Fig. 3 is the surface tension scaled by the factor IjQ, numerically evaluated using 
Eq. (29), as a function of —(1 — do)/Q. It is a monotonically decreasing function which goes to 
zero at the critical point. The agreement with the Laplace law was verified by simulation on a 
128 X 128 lattice. For initial conditions, we set up a circular high-density region in the center of the 
domain. When the system reaches equilibrium, a liquid bubble in very good circular shape forms 
and the bubble radius and the pressure difference inside and outside the bubble are then measured. 
According to the Laplace law, for a 2-D circular bubble, the pressure difference pi^ — Pout and the 
surface tension a are related by 

a 

Pin ~ Pout = ~^, (32) 

where R is the radius of the bubble. Plotted in Fig. 4 are the measurements and straight lines drawn 
withtheslopesgivenby the numerical solution of Eq. (29). The two cases are for —(l — c?o)/^ = 1.3 
and 1.4 respectively. The agreement between the simulation and the analytical results is good. 

As a simple application, we compute the dispersion relation of the capillary waves. In the 
absence of gravity, at long wavelength limit, the dispersion relation of capillary wave on a free 
surface is (see, e.g., Ref. [22]): 

co^ = (^] k\ (33) 

where a is the surface tension and p is the density of the fluid. The dispersion relation given by 
the lattice Boltzmann model was measured by examining a standing wave at a liquid-gas interface. 
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The simulations were carried out on a rectangular domain, measuring L x H in lattice unit, with 
periodic boundary conditions in both directions. Initially a interface parallel to the shorter edge 
was set up in the middle of the domain. A single sinusoidal wave with wavelength L was imposed 
on the interface and the amplitude of it was subsequently measured. To make the effects of finite 
water depth negligible, the aspect ratio of the domain, H/L has to be large and was chosen to be 
4 : ^ here. The maximum value of the wave number k = 2ir/L that can be simulated is limited by 
the fact that the wave amplitude has to be much larger than the thickness of the interface but much 
smaller than L. Four simulations with wavelength L = 256v^, 128-\/3, 64-\/3 and 32-\/3 were 
performed. The angular frequencies, cu, are plotted in Fig. 5 versus wave numbers k, together with 
Eq. (33). From the plot, the comparison is quite satisfactory. In these simulations, the parameters 
are — (1 — c?o)/^ = 1.1, c?o = 0.5, and r = 0.6. The surface tension and the densities of the two 
phases are a = 0.40, rig = 0.063, and ni = 2.23. 

V. DISCUSSION 

We have described in detail a lattice Boltzmann model for simulating fluids obeying a non-ideal 
gas equation of state. The equation of state can be changed arbitrarily. With a properly chosen 
equation of state, the LBE model can undergo a liquid-gas phase transition. In addition to the 
equation of state, we gave the bulk densities and the density profile at equilibrium, all as functions 
of a single temperature-like parameter. The surface tension is also given by numerical integration 
and it can be varied independent of the above equilibrium properties. 

This LBE model provides an efficient method for simulating flows involving interfaces and 
phase transitions. Phenomena near a liquid-gas critical point can also be simulated. The advantage 
over other LGA or LBE models developed for similar purpose are apparent, since in addition to 
its computational efficiency, almost everything with this model is known and can be changed to fit 
any desired properties. 

The collision time, r, should be made to depend on the local density when simulating a fluid 
system with large density variation, so that the two phases will have different transport coefficients. 
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The equilibrium properties discussed in this paper will not be altered by this change. 

As we have pointed out before, the long-range interactions in systems with multiple components 
can also be treated in LBE model using the same scheme. The phase diagrams of such a system 
are more complicated and remain to be calculated. In general, work similar to that contained in 
this paper can be carried out for the multi-component system as well. This is necessary before this 
method can be applied to practical problems to produce accurate quantitative results. 

The major weak point of this model, we believe, is the lack of an energy conservation relation 
and a dynamic temperature equation, in spite of the fact that a static temperature can be identified 
with a particular choice of the effective mass. Although lattice Boltzmann models with more 
speeds can be constructed to enable energy conservation [21], it is difficult to incorporate into it, 
at least with the scheme discussed here, the correct non-local inter-particle forces and conserve the 
total energy at the same time. This is mostly due to discrete lattice effect, that each time step is 
separated into a "streaming" and a "collision" step. It is difficult to keep the total energy conserved 
in these separated steps. 
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FIGURES 
FIG. 1. Coexistence curve when the function i/; is chosen to be 1 - exp{-n). The critical point is at 

— (1 — do)/Q = l.Sandric = ln2 ~ 0.693. The solid line is the solution of Eq. (25) and the dashed line is 

the solution to Eq. (30), the Maxwell construction. Both of them are obtained numerically. The diamonds 

are the results of numerical simulations. 

FIG. 2. The density profiles across the liquid-gas interface. Solid lines are solved from Eq. (24) by 
numerical integration. For each value of —(1 — do)/Q, results of two simulations with the angle between 
the normal vector of the interface and the lattice links being 30° and 0° are shown. No significant difference 
is found between them. As —(1 — do)/Q approaches the critical value, the distinction between the two 
phases becomes negligible. 

FIG. 3. Surface tension,scaledby 1/C/, as functionof — (l — c?o)/^- This curve is obtained by numerical 
integration using Eq. (29). 

FIG. 4. Laplace law tested by numerical simulations. The marks indicate measured pressure differences 
inside and outside a 2-D circular bubble pin — Pout , versus the reciprocal of the measured bubble radii. The 
lines are the analytic predictions. The two cases are for —(1 — do)/Q = 1.3 and —(1 — dQ)/Q = 1.4. 

FIG. 5. Dispersion relation of capillary wave at a liquid-gas interface. The diamonds are measured 
from lattice Boltzmann simulations and the dashed straight line is given by Eq. (33). 
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